intra-subject persistence for bacteria and virus SGBs

load packages

library(ggplot2)
library(data.table)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
data.table 1.14.2 using 1 threads (see ?getDTthreads).  Latest news: r-datatable.com
**********
This installation of data.table has not detected OpenMP support. It should still work but in single-threaded mode.
This is a Mac. Please read https://mac.r-project.org/openmp/. Please engage with Apple and ask them for support. Check r-datatable.com for updates, and our Mac instructions here: https://github.com/Rdatatable/data.table/wiki/Installation. After several years of many reports of installation problems on Mac, it's time to gingerly point out that there have been no similar problems on Windows or Linux.
**********

library(stringr)
library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:data.table’:

    between, first, last

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(tibble)
library(here)

here() starts at /Users/u241374/mike_tisza/teddy_virome/teddy_vir_bac_marker_geneR

library(rprojroot)
library(ggpubr)
library(rstatix)

Attaching package: ‘rstatix’

The following object is masked from ‘package:stats’:

    filter

set paths and filenames

### files
long_table=sprintf("%s/data/mp142_TGVG1.1_MPA4_combined_abundance_table_longform1.tsv", find_rstudio_root_file())
metadata_table=sprintf("%s/data/some_teddy_MP142_metadata2.all_samples1.delivery.csv", find_rstudio_root_file())

load long table and metadata, merge

long_dt <- fread(sprintf("%s", long_table), sep = "\t", header = T) %>%
  select(sampleID, rel_abundance, lineage) %>%
  mutate(kingdom = case_when(grepl("k__Bac", lineage) ~ "Bacteria", 
                             grepl("k__Vir", lineage) ~ "Virus",
                             grepl("k__Ar", lineage) ~ "Archea",
                             grepl("k__Euk", lineage) ~ "Eukaryota",
                             TRUE ~ "other"))

meta_dt <- fread(sprintf("%s", metadata_table), sep = ",", header = T) %>%
  select(-V1)

mask_id_count <- meta_dt %>%
  group_by(mask_id) %>%
  filter(n() >= 10) %>%
  summarize(total_samples = n())

merge_dt <- merge(long_dt, meta_dt, by.x = "sampleID", by.y ="sample")

persist_dt <- merge_dt %>%
  group_by(lineage, kingdom, mask_id) %>%
  summarize(samples_detected = n())

`summarise()` has grouped output by 'lineage', 'kingdom'. You can override using the `.groups` argument.

persist_dt <- merge(persist_dt, mask_id_count, by = "mask_id") %>%
  mutate(frequency_detected = samples_detected/total_samples)

dim(
  persist_dt %>% distinct(mask_id)
)

[1] 555   1

calculate stats

stat.test <- persist_dt %>% 
  filter(kingdom == "Virus" | kingdom == "Bacteria") %>%
  wilcox_test(frequency_detected ~ kingdom) %>%
  add_significance(p.col="p", output.col="p.signif")
stat.test <- stat.test %>% add_xy_position(x = "kingdom")
stat.test

NA
NA

plot violins

persist_dt %>% 
  filter(kingdom == "Virus" | kingdom == "Bacteria") %>%
  ggplot(aes(factor(kingdom), frequency_detected)) + 
  geom_violin(aes(fill = factor(kingdom)),  alpha = 0.6, draw_quantiles = c(0.25, 0.75), linetype = "dotted") +
  geom_violin(fill="transparent", draw_quantiles = c(0.5)) +
  theme_bw() + 
  xlab("") + 
  ylab("SGB Intra-subject Frequency") + 
  scale_fill_manual(values=c("#F2AD00", "#5BBCD6")) + 
  stat_pvalue_manual(stat.test, label = "p.signif", y.position = 1.05, tip.length = 0.01, coord.flip = TRUE)  + 
  coord_flip() +
  theme(legend.position = "none")

ggsave(sprintf("%s/charts/bacteria_vs_virus_SGB_intrasubject_frequency1.pdf", find_rstudio_root_file()), width = 4, height = 3)

stat.test <- as.data.frame(stat.test)
stat.test <- apply(stat.test,2,as.character)
write.csv(stat.test, file=sprintf("%s/intermediate_files/wilcoxon_test_vOTU_bOTU_persist_10samps.csv", find_rstudio_root_file()))